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Abstract 

In practice, several time series exhibit long-range dependence or per- 
sistence in their observations, leading to the development of a number of 
estimation and prediction methodologies to account for the slowly decaying 
autocorrelations. The autoregressive fractionally integrated moving average 
(ARFIMA) process is one of the best-known classes of long-memory models. 
In the package afmtools for R, we have implemented some of these statistical 
tools for analyzing ARFIMA models. In particular, this package contains 
functions for parameter estimation, exact autocovariance calculation, predic- 
tive ability testing, and impulse response function, amongst others. Finally, 
the implemented methods are illustrated with applications to real-life time 
series. 

Key words: ARFIMA models, long-memory time series, Whittle esti- 
mation, exact variance matrix, impulse response functions, forecasting, R 
package. 



1 Introduction 



Long-memory processes introduced by Granger and Joyeux ( 1980[ ) and Hosking 
(1981), are playing a key role in the time scries series literature (for example 



see Palma (2007) and references therein) and have become a useful model for 
explaining natural events studied in geophysics, biology, and other areas. As a 
consequence, a number of techniques for analyzing these processes have been de- 
veloped and implemented in statistical packages. For example, packages about 



long- memory processes have been developed in R (R Development Core Team 



2012): the longmemo package produces a Whittle estimation for Fractional Gaus- 



sian Noise and Fractional APJMA models via an approximate MLE using the 



Beran (19941 algorithm and performs Spectral Density of Fractional Gaussian 
Noise and Periodogram Estimate. In addition, the fracdiff package simulates 
ARFIMA time series, estimates ARFIMA parameters using an approximate MLE 
approach (Haslett and Raftery 19891, and calculates their variances with the Hes- 
sian method. Recently, |Hyndman and Khandakar (20081 describe the forecast 
package to automatically predict univariate time series via State Space models 
with exponential smoothing for ARIMA models. In addition, the forecast package 
offers a forecast function for ARFIMA models estimated using the algorithm pro- 
posed by |Peiris and Perera (1988). The afmtools package requires the polynom, 
hypergeo, sandwich and the aforementioned fracdiff and longmemo packages. 

Unfortunately, many of these computational implementations have important 
shortcomings. For instance, there is a severe lack of algorithms for calculating ex- 
act autocovariance functions (ACVF) of ARFIMA models, for computing precise 



estimator variances, and for forecasting performance tests (Giacomini and White 



2006), and impulse response functions ([Hassler and Kokoszka 2010), as well as 



for other aspects. In order to circumvent some of these problems, this paper dis- 
cusses the package afmtools developed by |Contreras-Reyes et ah (2011). This 
package aims to provide functions for computing ACVFs by means of the |Sowell| 
( 1992 1 algorithm, ARFIMA fitting through an approximate estimation scheme via 
Whittle algorithm (Whittle 19511, asymptotic parameter estimate variances and 
several other tasks mentioned before. Hence, the aims of this paper are to analyze 
the afmtools package and to illustrate its theoretical and practical performance, 
which complements the existing development packages related to ARFIMA mod- 
els mentioned above. Specifically, we implement our findings in a meteorological 
application about tree ring growth. 

The remainder of this paper is structured as follows. Section 2 is devoted 
to describing the ARFIMA processes and their properties. This section includes 
an analysis of the spectral density, autocovariance function, parameter variance- 
covariance matrix estimation, impulse response function, and a model parameters 
estimation method. In addition, this section provides a test for assessing the pre- 
dictive ability of a time series model. Finally, Section 3 addresses the performance 
of the functions implemented in the afmtools package. Apart from describing the 
methodologies implemented in this package, we also illustrate their applications 
to real-life time series data. 



2 ARFIMA processes 

Recent statistical literature has been concerned with the study of long-memory 
models that go beyond the presence of random walks and unit roots in the uni- 
variate time series processes. The autoregressive fractionally integrated moving- 



average (ARFIMA) process is a class of long-memory models ( Granger and Joyeux 
p80| ; |Hosking| ( p8lj )), the main objective of which is to explicitly account for 



persistence to incorporate the long-term correlations in the data. The general 
expression for ARFIMA processes {yt} may be defined by the equation 

®(B)y t =e(B)(l-B)- d e t , (1) 

where $(5) = 1 - faB (f> p BP and 6(5) = 1 + 9 ± B H V B q Bi are the 

autoregressive and moving-average operators, respectively; Q(B) and <d(B) have 
no common roots, B is the backward shift operator and (1 — B)~ d is the fractional 
differencing operator given by the binomial expansion 

(1 - B)- d = Y TU ± d \ , B* = Y m& , (2) 
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for d € (—1,1/2) and {et} is a white noise sequence with zero mean and 
innovation variance a 2 . An asymptotic approximation of 



for large j is 



„ - r ^ + d ) ( 3 ) 



jd-l 

* ~ wr (4) 



where T is the usual gamma function. 

Theorem 2.1 Consider the ARFIMA process defined by (JT|) and assume that the 
polynomials $(■) and &(■) have no common zeros and that d £ (—1, |)- Then, 

a) // i/ie zeros o/ <&(•) Zze outside the unit circle {z : |z| = 1}, i/ien iftere is a 

unique stationary solution of (JlJ given by yt — X)^L-oo ' l Pj £ t-j where tpj are 
the coefficients of the following polynomial tp(z) = (1 — z)~ d Q(z)/$>(z). 

b) If the zeros of <&(•) lie outside the closed unit disk {z : \z\ < 1}, then the 

solution {y t } is causal. 

c) If the zeros of &(■) lie outside the closed unit disk {z : \z\ < 1}, then the 

solution {yt} is invertible. 



For a proof of Theorem 2.1 see e.g. Palma (2007). Recall that, according to 
the representation theorem of Wold (19381, any stationary process is the sum of 
a regular process and a singular process; these two processes are orthogonal and 
the decomposition is unique. Thus, a stationary purely nondetcrministic process 
may be expressed as 

oo 

y t =ip(B)et=^ri> j £t-j, (5) 

3=0 

The spectral measure of the purely nondeterministic process ^ is absolutely 
continuous with respect to the Lebesgue measure on [— tt,tt], where the spectral 
density of the process ([T]) can be written as 



/(A) = r l#- 

ZTT 



(G) 



a 
2tt 
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where i denotes the imaginary unit. A special case of ARFIMA models is 
ced pr< 

Q(B) — 1 and the spectral density is given by 

/(A) = 



the fractionally differenced process described by Hosking (1981), in which the 
polynomials are <&(£?) 



-i\\—2d 



2tt' 
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2.1 Whittle estimation 



The methodology to approximate MLE is based on the calculation of the pe- 
riodogram I(X) by means of the fast Fourier transform (FFT); e.g., Singleton 
( |1979| ), and the use of the approximation of the Gaussian log-likelihood function 
due to Whittle (19511 and by Bisaglia and Guegan (19981. So, suppose that the 



sample vector Y = (yi, y2, ■ ■ ■ , y n ) is normally distributed with zero mean and 
autocovariance given by ( 16 1 as 



j(k - j) 



f{\)e lX{k - ]] d\ 



(7) 



where /(A) is defined as in ((6| and is associated with the parameter set SI of 
the ARFIMA model defined in (HI. The log likelihood function of the process Y 
is given by 



i(^) = -^[log|A|-Y T A- 1 Y] 



(8) 



where A = [y(k — j)] with k,j = 1, ...,n. For calculating pj, two asymptotic 
approximations are made for the terms log(|A|) and Y T A _1 Y to obtain 



L(U) w - 



1 

-in 



log[27r/(A)]dA 



dX 



(9) 



as n — y oo, where 



/(A) 



1 

2im 



(10) 



is the periodogram indicated before. Thus, a discrete version of is actually 
the Riemann approximation of the integral and is 



L(fl) w - 



1 

2n 



^io g /(A,)+x: 



/(A,-) 



(11) 



where Xj — 2i:j/n are the Fourier frequencies. Now, to find the estimator 
of the parameter vector fi, we use the minimization of L(f2) produced by the 
nlm function. This non-linear minimization function carries out a minimization 
of L(Tl) usi ng a Newton-type algorithm. Under regularity conditions according 
to Theorem 2.2 (see Section 2.2), the Whittle estimator SI that maximizes the 



log- likelihood function given in (11) is consistent and distributed normally (e.g. 
|Dahlhaus[[l994t . The following fi gures illustrates the performance of the Whittle 
estimator for an ARFIMA(1, d, 1) model. 

Figure [l] shows several simulation results assessing the performance of the es- 
timators of d, AR, and the MA parameters, for different ARFIMA models. These 
plots include the exact and Hessian standard deviations. According to the defini- 
tion of the ARFIMA model, the simulations are run in the interval (—1, 0.5) for d. 
In addition, Figure [2] shows some simulation results regarding the log-likelihood 
behavior for the cases d = {—0.9, —0.6, —0.3, 0, 0.25, 0.45} with a rectangular grid 
x 9 = (-0.9,0.9) x (-0.9,0.9) for the ARFIMA ( 1, d, 1) model. The plots of 
Figure [T] show a similar behavior for the estimators with respect to the theoretical 
parameters, except for the extreme values of the ARMA parameters near -1 and 
1. Consequently, the confidence intervals tend to be larger than the other values 
of 4> an d parameters for plots (b) and (e) . In addition, the plots of Figure [2] 
present low values of the likelihood function for the values of (f> and 9 closed to 0, 
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Figure 1: Dispersion plots between estimated and theoretical parameters of (a): 
ARFIMA(0,d,0), (b)-(c): ARFIMA(1, d, 0) and (d)-(e): ARFIMA(0, d, 1) where 
the green dotted line is the exact standard deviation and the red line is the Hessian 
standard deviation. 



especially for the plots (c)-(f) when d = {—0.3,0,0.25,0.45}. However, the plots 
(a)-(d) show high values of the likelihood function when this is evaluated for the 
points near (j> — —0.9 and 9 — 0.9. For the plots (e)-(f), the behavior is inverse, i.e., 
the likelihood function tends to be higher for values near <j> = 0.9 and 8 = —0.9. 



2.2 Parameter variance-covariance matrix 

Here, we discuss a method for calculating the exact asymptotic variance-covariance 
matrix of the parameter estimates. This is a useful tool for making statistical 
inferences about exact and approximate maximum likelihood estimators, such as 



the Haslett and Raftery ( 1989 1 and Whittle methods (see Section 2.1). An example 
of this calculation for an ARFIMA(1, d, 1) model is given by Palma (2007 pp. 105- 
108). This calculation method of the Fisher information matrix is an alternative 
to the numerical computation using the Hessian matrix. 

This proposed method is based on the explicit formula obtained by means of the 
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Figure 2: Log-Likelihood 3D plots for ARFIMA(1, d, 1) model with parameters 
(a) d = -0.9, (b) d = -0.6, (c) d = -0.3, (d) d=0,(e)d= 0.25 and (f) d = 0.45 
using a grid 0X0 = (-0.9, 0.9) x (-0.9, 0.9). 

derivatives of the parameters log-likelihood gradients. From the spectral density 
defined in we define the partial derivatives V$ = (d/dfa) and Ve = (d/d0j), 
with i = 1, . . . , p and j = 1, . . . , q. 

Theorem 2.2 Under the assumptions that y t is a stationary Gaussian sequence, 
the densities f (A) , / _1 (A), d/d^if~ 1 {\), d 2 /dfiid[ijf~ 1 (A) and d 3 /d^id^jd^kf^ 1 (^) 
are continuous in (A, fx) for a parameter set /i = {d, 4>i, ■ ■ ■ , 4>p, 0%, ■ ■ ■ , q }; we have 
the convergence in distribution for an estimated parameter Jl and the true param- 
eter Ho about a Gaussian ARFIMA model with 

Vnfo-iio) ^(O,^- 1 ^)), (12) 



G 



where 



= / [V%/ M (A)][V%/ M (A)] T dA. 



(13) 



For a proof of Theorem 2.2 see e.g. Palma (20071. Thus, if we consider the 
model with spectral density ^ where {et} is an independent and identically 
distributed iV(0,<7 2 ), we have that the parameter variance-covariance matrix T 
may be calculated in the following proposition. 



Proposition 2.3 If {y t } is stationary, then 
d 

m 
d 



logf(X) = -%[2(l-cosA)], 

E?=i^cos[(^-j)A] 



-logf(X) 



d 



logf(X) = 



E?=i ELi 4>ifa cos i(j - fc ) A l 

E?=i^cos[^-i)A] 



Mi ^ J v E?=i ELi cos[(j - A)A] ' 
Proof First, from the spectral density given in ([6| we have that 



log/(A) = log ( — ) - <flog[2(l - cos A)] + log|6( e jA )| 2 - log|0>(e 



\\ |2 



By Theorems 



2.1 



and 



2.2 



iA\|2 



we observe that $(e lA ) = Ej=i ^j^K this yields 

j=i fe=i 
p P 

= 2££^0 fc cos[O--fc)A], 
j=i fc=i 



and 



|$(e 



i\\\2 



i\{l-k) 



P V 



j=i k=i 
p 



3=1 



Analogously, we have that 



30, 



|e( e iA )| 2 = 2£^co S [(^-i)A]. 

3=1 



Then, this implies the results for g^jlog/(A) and -^log/(A). For ^log/(A) is 
direct . 

Some computations and implementation of this matrix are described in Sec- 
tion 3.2, associated with the parameters of several ARFIMA models, using the 
Whittle estimator. 
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2.3 Impulse response functions 



The impulse response functions (IRF) is the most commonly used tool to eval- 
uate the effect of shocks on time series. Among the several approximations to 
compute this, we consider the theory proposed by Hassler and Kokoszka| (20101 
to find the IRF of a process {y t } following an ARFlMA(p,d,q) model. The 
properties of these approximations, depend on whether the series are assumed 
to be stationary according to Theorem |2.1| So, under the assumption that the 
roots of the polynomials Q(B) and <d(B) are outside the closed unit disk and 
d S (—1, 1/2), the process {y t } is stationary, causal and invertible. In this case, 
we can write y t = ^{B)e t where ^f(B) represents the expansion of the MA(oo) 
coefficients denoted as ipj with j > 1. These coeffi cients satisfy the asymptoti c 



relationship ipj ~ [e(l)j' d " 1 ]/[$(l)r(d)] as j -4 oo ( |Kokoszka and Taqqu[ [l995[ ), 
6(1) = 1+X}i=i ^ij an d $(1) = 1 — Ef=i As a particular case, we have that the 
ipj coefficients for an ARFIMA(0, d, 0) are given in closed form by the expression 
ipj = T(j + d)/(T(j + l)r(d)). Now, from (pi and the Wold expansion the 



process (jl| has the expansion (1 — B) y t = z2jLo r ljyt-j — J^jLo Rj £ t-j> where 
Rj is the so-called IRF and is given by 

j 

i=0 

The terms rjj can be represented in recursive form using ([3])asm = (1 + (d — l)/j)r)j—i, 
for j > 1 and ijq = 1. From the asymptotic expression given in Q and assuming 
that Y^jLo < °°i we have the following asymptotic representation 

3 d - 1 ^ 



as j — > oo and i/jj/(j d — ► 0. 



2.4 Autocovariance function 



We illustrate a method to compute the exact autocovariance function for the gen- 
eral ARFIMA(p, d, q) process. Considering the parameterization of the autocovari- 
ance function derived by writing the spectral density ([6| in terms of parameters 
of the model given by Sowell (19921, the autocovariance function of a general 



ARFIMA(p, d, q) process is given by 



7(A) 



1 

2^ 



/(A)e 



-iXh 



dX, 



(16) 



where i denotes the imaginary unit. Particularly, the autocovariance and au- 
tocorrelation functions of the fractionally differenced ARFIMA(0, d, 0) process are 
given by 



70 (h) 
Po(h) 



2 r(l-2d) T(h + d) 
° r(l - d)T(d) T(l + h-d)' 
r(l-d) T(h + d) 
T(d) T(l + h-d)' 
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respectively. Then, the polynomial Q(B) in (|T|) may be written as 

v 

$(B) = H(l- Pi B). (17) 

Under the assumption that all the roots of 4>(B) have multiplicity one, it can 
be deduced from (16) that 



q p 



i=—q j=l 



with 



Hi = 



Pi 



- piPj)Y[(Pj - pk) 



C(d,h,p) 

13(h) 

F(a, 6, c, x) 



-^[p*P m + P{-h)-l], 



F(d + h, 1, 1 - d + h,p), 

a ■ (a + 1) • b ■ (b+ 1) 2 
c-(c+l)-l-2 X 



a ■ b 

1 + — -x 
c ■ 1 



where F(a, b, c, x) is the Gaussian hypergeometric function (e.g. Gradshteyn 
and Ryzhik{|2007[ ). The term pr esented here a nd in |Palma| ( |2007[ pp. 47-48) 
is a corrected version with respect to Sowell ( 1992 ) and is 



min((7,g+z) 

ip(i) = OkOk-i- 

k— max(0,i) 

In the absence of AR parameters the formula for j(h) reduces to 

a « r(i-2d)r(fe + d-i) 



On the other hand, the findings of Hassler and Kokoszka (2010 1 describe the 



asymptotic behavior of the autocovariancc function "f(h) as 



7 (/i) ~ Cj \h\ 2d -\ 



(18) 



where c 7 = a 2 n l r(l-2d) sin(-7rd) fejlo ^j) for lar S e N- Let {z/i, y 2 , - • ■ , J/n} 
be a sample from the process in (|T|) and let y be the sample mean. The exact vari- 
ance of y is given by 



Var(y) 



J 



n— 1 / 



By (18 1 and for large n, we have the asymptotic variance formula Var(y) ~ 
n 2d ~ 1 c~ f /d(2d+l). The Sowell method is implemented in Section 3.3 for a selected 
ARFIMA(1, d, 1) model and its exact autocorrelation function is compared with 
the sample autocorrelation. 
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2.5 Predictive ability test 



One approach to compare prediction models is through their root mean square 
error (RMSE). Under this paradigm and given two forecasting methods, the one 
that presents the lower RMSE is the better of the two. To compare statistically 
the differences of predictive ability among two proposed models, we focus here on 
the evaluation paradigm proposed by Giacomini and White (2006. GW). This test 
aims to evaluate a prediction method but not to carry out a diagnostic analysis. 
Therefore, it does not consider the parametric uncertainty, which is useful if we 
want to compare nested models from an ARFIMA model. The GW test attributed 
Diebold and Mariano ( 1995 1 is based on the differences AL, ; = \xi 



to 



-Vi 



-y, 



where Xi and % are the forecasted observations of the first and second model 
respectively, for i — 1, n. The null hypothesis for the GW test associated with 
expected difference _E[AL] of the two prediction models is Ho'- E[AL] = 0, whereas 
the alternative is H\\ E[AL] ^ 0. These hypotheses are tested by means of the 
statistic, 

AL(N) -- ' 




A Li 



where N — 



to + 



1, n is the total size of the sample, r is the prediction 
horizon, and to is the observation at which the mobile windows start. Note that 
under Ho, the statistic AL(N) is asymptotically normal. For r — 1, an estimator 
of (7/v can be obtained from the estimation of a-^ from a simple regression AL(N) — 
2 + e. However, for horizons r > 1, it is possible to apply a heteroscedasticity and 
autocorrelation consistent (HAC) estimator; for example, Newey and West ( |1987 1 



Andrews (1974). The GW test is implemented in the gw.testO function and 



or 



described in Section 3.4. 



3 Application 

The functions described in the previous sections arc implemented in the af mtools 
package. We illustrate the performance of the Whittle method by applications to 
real-life time series TreeRing (Statlib Data Base, http://lib.stat.cmu.edu/) dis- 
played in Figure [3] left. In climatological areas, it is very important to analyze this 
kind of data because this allows us to explore rainy and dry seasons in the study 
area. Hipel and McLeod (19941 has been modeling this time series to determine 
the range of possible growths for the upcoming years of the trees using ARMA and 
ARIMA models. On the other hand, this time series displays a high persistence in 



its observations and has been analyzed by Palma and Olea (20101) and Palma et 



al. (20111 with a locally stationary approach. Apparently, the illustrated growth 
of the trees, represented by the number of the rings, presents a long-range depen- 
dence of its observations along the observations for ages and seasons (see Figure [3j 
right). For these reasons, we model the Tree Ring widths time series using long- 
memory models; specifically, the ARFIMA models are used to estimate, diagnose, 
and compare forecasts of the number of tree rings for upcoming years. 

In order to illustrate the usage of package functions, we consider a fitted 
ARFIMA(1, d, 1) model. For this model, we have implemented the Whittle al- 
gorithm and computed the exact variance-covariance matrix to compare with the 
Hessian method. Afterward, we compare the Sowell method for computing the 
ACVF function with the sample ACVF. Other functions have also been imple- 
mented and illustrated in this section. 
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Figure 3: Left: Tree Ring Data Base. Right: Illustration of the observed climatic 
episodes and tree ring growth. Image source: http://summitcountyvoice.com 



3.1 The implementation of the Whittle algorithm 

The estimation of the fractionally, autoregrcssivc, and moving-average parameters 
has been studied by several authors (Beran (1994); Haslett and Raftery] ( |1989[ ); 
Hyndman and Khandakar (2008)). A widely used method is the approximate 
MLE method of iHaslett and Raftery! (119891). In our study, estimation of the 



ARFlMA(p,d,q) model using the corresponding Whittle method (Whittle 1951) 



is described in Section 2.1 and this model is fitted by using the arfima. whittle 
function. To apply the Whittle algorithm to the TreeRing time series as an exam- 
ple, we use the following command considering an ARFIMA(1, d, 1) model. 



R> y <- data (TreeRing) 

R> model <- arfima. whittle (series = y, nar = 1, nma = 1, fixed = NA) 

Note that the option fixed (for fixing parameters to a constant value) has 
been implemented. This option allows the user to fix the parameters d, cf>i, or 6\, 
in order of occurrence. For example, in our ARFIMA(l,d, 1) model, we can set 
the parameter d to be equal to zero. Consequently, we obtain the estimation of a 
simple ARMA(1, 1) model. The object model is of class arfima and provides the 
following features: 

• estimation of d and ARM A parameters; 

• standard deviation errors obtained by the Hessian method and the respective 
t value and Pr (> 1 1 1 ) terms; 

• the log-likelihood function performed in the arfima. whittle . loglikO 

function; 

• innovation standard deviation estimated by the ratio of the theoretical spec- 
trum; 

• residuals from the fitted model. 

The commands plotO, residualsO, summaryO, tsdiagO and printO 

have been adapted to this model class of S3 method. The summaryO option 
shows the estimated parameters, the Hessian standard deviations, the i-statistic, 
and their respectively p- values. The computation of the long- memory parameter 
d as well as the autoregressive {4>\, <p p } and moving average {0\, ■■■,0 q } param- 
eters can be handled quickly for moderate sample sizes. Printing the model object 
by the summaryO function shows the items mentioned before as 
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Table 1: Summary of estimated parameters for several ARFIMA models. 







FN(d) 




ARFIMA(1, 


i,0) 


AR.FIMA(0, d, 1) 


ARFIMA(1, d 


,1) 


Par. 


Est 


Hessian 


Exact 


Est Hessian 


Exact 


Est Hessian Exact 


Est Hessian 


Exact 


d 


0.195 


0.048 


0.023 


0.146 0.048 


0.038 


0.156 0.048 0.035 


0.106 0.048 


0.063 
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0.072 0.029 


0.049 


0.059 0.029 0.045 


0.397 0.035 
-0.285 0.032 


0.282 
0.254 



R> summary (model) 
$call 

arf ima. whittle (series 



1) 



$coef mat 

Estimate 
d 0.1058021 
phi 1 0.3965915 
theta 1 -0.2848590 



Std. Error t value Pr(>|t|) 
0.04813552 2.198004 0.02794879 
0.03477914 11.403142 0.00000000 
0.03189745 -8.930462 0.00000000 



$sd. innov 
[1] 35.07299 

$method 

[1] "Whittle" 

attr( , "class") 

[1] "summary . arf ima" 

Furthermore, we evaluate the Whittle estimation method for long- memory 
models by using the theoretical parameters versus the estimated parameters (see 
Figure [TJ. In addition, we compare the exact standards deviations of the param- 
eters from Section 2.2 with those obtained by the Hessian method. These results 
and the evaluation of the Whittle estimations are illustrated in Table [T] 



3.2 The implementation of exact variance-covariance ma- 
trix 

The var . af m() function shows the exact variance-covariance matrix and the stan- 



dard deviations. The computation of the integrals of the expression (13 1 is carried 



out by using the Quadpack numeric integration method (Piessens et al. 1983) 



implemented in the integrate () function (stats package). Note that the func- 
tions involved in these integrals diverge in the interval A = [— 7T, 7r]. However, 
they are even functions with respect to A. Thus, we integrate over [0,7r] and then 
multiply the result by two. Now, by using the central limit theorem discussed 
in Section 2.2, we can obtain the asymptotic approximation of the estimated pa- 
rameters standards errors SE(«)i = (ip" 1 ],,) 1 / 2 of an ARFIMA model, where 
f2 = (d, 0i, . . . , tfip, 0i, ... , 6 q ) and [E^ 1 ]^ corresponds to the zth diagonal compo- 
nents of the matrix S _1 for i — {1, ...,p + q + 1}. 

By using the Whittle estimators, we search for the lowest AIC (Akaikc Infor- 
mation Criterion, Akaike 1974) given by AIC(w) = — 2[logi(S) — (p + q+ 1)] over 
a class of ARFIMA models with p, q G {0, 1, 2}, where w is a subset of f2 and L(uj) 
is the likelihood associated with w. From Table [2] we can see that the fraction- 
ally differenced model ARFIMA(0, d, 0) has the lowest AIC. Candidate models are 
marked in bold in Table [21 
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Table 2: Akaike's criterion for £1 = (d, <fii, <p2, 01,82) with p- values obtained for the 
Hessian standard deviation. 



p 


q 


a ir* 

AILj 


a 


p- value 








-37.018 


0.196 








1 


-35.008 


0.156 


0.001 





2 


-33.000 


0.113 


0.018 


1 





-35.007 


0.146 


0.002 


1 


1 


-33.004 


0.106 


0.028 


1 


2 


-30.996 


0.142 


0.003 


2 





-33.002 


0.111 


0.021 


2 


1 


-30.995 


0.130 


0.007 


2 


2 


-28.915 


0.191 






Standardized Residuals 




ACF of Residuals 



p values for Ljung— Box statistic 
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Figure 4: Plots of residuals analysis from tsdiag standard command adapted to 
ARFIMA model residuals. 



Additionally, we propose a technique for obtaining the spectral density associ- 
ated with ARFIMA and ARMA processes in spe ctrum . ar f ima ( ) and spe ctrum . arma ( ) , 
respectively. This is done by using the polyrootO function of the polynom pack- 
age to compute the roots of the polynomials <!>(e -lA ) and 0(e _iA ). Both functions 
need the estimated ARFIMA parameters and the sd.innov innovation standard 
estimation given by an object of arf ima class. For the spectrum density and pe- 
riodogram, see Section 2 and Subsection 2.1, respectively. Since the calculation 
of the FFT has a numerical complexity of the order 0[nlog 2 (n)], this approach 
produces a very fast algorithm to estimate the parameters. It is possible to obtain 
the FFT through the fft() function based on the method proposed by Singleton 
(II979I. 
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Inverse AR and MA Roots 



ACF 




Figure 5: Diagnostic Plots made by the plot () command. Left: Plots of unit roots 
circle along with the root provided by AR and MA polynomials and Theoretical 
(red line) vs. Empirical Spectrum (back points) plot. Right: ACF plots of Tree 
Rings and ARFIMA(1, d, 1) model residuals. 



3.3 The implementation of the diagnostic functions 

We have also implemented a very practical function called check . parameters . arf ima 
This verifies whether the long-memory parameter d belongs to the interval ( — 1, 0.5) 
and whether the roots of the fitted ARM A parameters lie outside the unit disk. 
This function was incorporated in the plotO command. In the first plot of Fig- 
ure [5j we can see that the roots of the AR and MA polynomials lie outside the 
unit disk, according to the assumptions of stationarity solutions of ([lj presented in 
Theorem2.1 (see Section 2). Alternatively, the check. parameters. arf ima() that 



takes an arf ima-class object, gives TRUE/FALSE- type results indicating whether the 
parameters pass the requirement for a stationary process. 

Additionally, an adaptation of the function tsdiagO can be found in this 
package. This is implemented in an S3- type arf ima class method and shows three 
plots for analyzing the behavior of the residual from the fitted ARFIMA model. 
This function has additional arguments such as the number of observations n for 
the standardized residuals and critical p- value alpha. Figure [4] illustrates these 
results, where, the residuals are white noise at a confidence level of a = 95%. 

On the other hand, the Rj (IRF) illustrated in Figure [7] decays exponentially 
fast, at a rate of because, these functions inherit the behavior of rjj. This 
behavior is typical for ARFIMA models, as reported by |Hassler and Kokoszka| 
( |2010| ), |Kokoszka and Taqqu| ( |1995[ ) and |Hosking| ( |1981| ). Figure [7j shows some Rj 
curves associated with the three models considered in Table [l] for the asymptotic 
method by formula (151 (labeled Asymptotic in the plot) and the counterpart 
method by formula ( 14 1 (labeled Normal in the plot). Note that for a large value 
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Figure 6: Theoretical (left) and Empirical (right) ACF of selected ARFIMA model. 
Impulse Response: h=50, ARFIMA(1.0.106.1) Impulse Response: h=50. ARFIMA(1.0.146.0) Impulse Response: h=50, ARFIMA(0,0.156,1) 
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Impulse Response: h=100. ARFIMA(1,0.106,1) Impulse Response: h=100. ARFIMA(1,0.146,0) Impulse Response: h=100, ARFIMA(0,0.156,1) 
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Impulse Response: h=150, ARFIMA(1.0.106.1) Impulse Response: h=1 50. ARFIMA(1,0.146,0) Impulse Response: h=150, ARFIMA(0,0.156,1) 
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Figure 7: Plots of IRFs methods for h = 50, 100 and 150 lags. 



of j » 50, both methods tend to converge, and, the curves make an inflexion in 
the value j ~ 10. Note that the Asymptotic approximation tends to be equal to 
the Normal method in the measure that the input h lag increases (see plots for 
h=150 in Figure 0. These IRFs are available in the function ir.arfimaO, with 
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arguments h to evaluate the IRFs over a particular h lag and, model for an object 
arfima. whittle. The ir.arfimaO function produces the vectors RE and RA for 
Normal and Asymptotic IRFs. 

The exact Sowell autocovariance computation obtained by rho.sowellO and 
the sample autocorrelation obtained by the ACF() command are applied to tree 
ring time series. In Figure [6] the blue dotted lines in the first plot correspond to the 
{=p2/y/n} significance level for the autocorrelations. The function rho.sowellO 
requires the specification of an object of class arfima in the object option that, by 
default, is NULL. But, if object=NULL, the user can incorporate the ARFIMA pa- 
rameters and the innovation standard deviation. Alternatively, the implemented 
plot option gives a graphical result similar to the ACF ( ) command in the sample 
autocorrelation. We can see the similarity of both results for the discussed model. 
The ACVF implementation is immediate but, for the calculation of the Gaussian 
hypergeometric functions, we use the hypergeoO function from the hypergeo 
package. For values of h > 50, we use the approximation (18 1 reducing consider- 
ably the computation time as compared to the Sowell algorithm. 

On the other hand, the rho.sowellO function is required by the function 
smv . af m ( ) . The function smv . af m ( ) calculates the variance of the sample mean of 
an ARFIMA process. When the argument comp is TRUE, the finite sample variance 
is calculated, and when comp is FALSE, the asymptotic variance is calculated. 




Figure 8: Plots of out-of-sample predictions for (a) r = 2, (b) r = 4 and (c) r = 5 
for two models: ARIMA(0,1,0) in red and ARFLMA(1, d, 1) in blue with Haslett 
& Raftery estimator (HR). 
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Table 3: Summary of p- values of the GW test for each HAC estimator, r = {2, 4, 5} 
prediction horizon parameters and for the estimator methods B (Benchmark 
model), ML (ARFIMA models using Maximum Likelihood estimator) and HR 
(ARFIMA models using Haslett & Raftery estimator) over 40 observations of sam- 
ples from the TreeRing data set. The p-values marked in bold are lower than the 
probability (0.05) related to a 5% confidence level. 
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0.018 
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HR 
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3.4 Forecasting evaluations 



The GW method implemented in gw.testO for evaluating forecasts proposed by 



Giacomini and White ( 2006 I compares two vectors of predictions, x and y, provided 
by two time series models and a data set p. We consider that it is relevant to 
implement this test to determine if the predictions produced by a time series model 
(e.g., ARFIMA) process good forecasting qualities. This test for predictive ability 
is of particular interest since it considers the tau prediction horizon parameter or 
ahead in the case of pred.arf ima() function. Alternative methods are discussed, 
for instance, by Diebold and Mariano (19951. If tau=l, the standard statistic 
simple regression estimator method is used. Otherwise, for values of tau larger 
than 1, the method chosen by the user is used in the method option. The available 
methods for selection are described below. They include several Matrix Covariance 
Estimation methods but, by default, the HAC estimator is used in the test. The 
user can select between the several estimators of the sandwich package mentioned 
before: 

• HAC: Heteroscedasticity and Autocorrelation Consistent (HAC) Covariance 
Matrix Estimation by vcovHACQ function ( |Zeileis| ( p2004| ); |Zeileis| ( |2006[ )) . 



NeweyWest: Newey-West HAC Covariance Matrix Estimation by NeweyWest () 



function (Newey and West (1987)). 



LumleyHeagerty: Weighted Empirical Adaptive Variance Estimation by 
weave () function (Lumley and Heagerty (1999])). 



• Andrews: Kernel-based HAC Covariance Matrix Estimation by kernHACO 



function (Andrews (1974); Andrews and Monahan (1992[)) 



This test gives the usual results of the general test implemented in such as the 
GW statistic in statistic, the alternative hypothesis in alternative, the p- value 
in p. value, others such as the method mentioned before in method, and the name 
of the data in data . name. In some studies, the GW test is used to compare selected 
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models versus benchmark models such as ARM A, ARIMA, or SARIMA models 



(e.g. Contreras-Reyes and Idrovo 2011). So, to illustrate the GW test perfor- 



mance, we simulate the out-of-sample prediction exercise through mobile windows 
for TreeRing data sets considering the first 1124 observations and, later, forecast- 
ing 40 observations using three models: the ARIMA(0,1,0) benchmark model, the 
ARFIMA(p, d, q) with MLE estimator, and the Haslett & Raftery estimator (HR) 
using the algorithm of the automatic forecast () function implemented in the 
forecast package by Hyndman and Khandakarj ( 2008 1 algorithm. In Figure [8j 
the GW test compares the out-of-sample predictions of the ARIMA(0,1,0) and 
ARFIMA(p,d,q) model. It is important to note that the goal of this test is only 
to compare prediction abilities between models. 

Finally, we study a more general simulation, comparing the three predictors 
vectors with 40 real observations using gw.testO function considering the hy- 
potheses testing alternative="two . sided" to contrast significant differences be- 
tween predictions. In addition, we consider the four HAC estimators mentioned in 
the beginning of this subsection and prediction horizon parameters r = {2,4,5}. 
The results are summarized in Table [3] and Figure [8] We can see that the differ- 
ences in the prediction ability between B vs. MLE and B vs. HR are significant 
for r = 2 and 4 but, between MLE and HR they are not unequal for the three 
considered values of r. Given the non-significance of the MLE-HR test p- value, 
the MLE model is not considered in Figure [8] As expected, the increments of the 
prediction horizon parameter showed that the different prediction abilities of the 
models tend to be zero because the time series uncertainty tends to increase. Con- 
sequently, the individual prediction performance of each model is not considered 
by the test. 



4 Conclusions 

We developed the af mtools package with the goal of incrementing the necessary 
utilities to analyze the ARFIMA models and, consequently, it is possible to exe- 
cute several useful commands already implemented in the long- memory packages. 
In addition, we have provided the theoretical results of the Whittle estimator, 
which were applied to the Tree Ring data base. Furthermore, we have performed 
a brief simulation study for evaluating the estimation method used herein and also 
have evaluated the properties of its log-likelihood function. The numerical exam- 
ples shown here illustrate the different capabilities and features of the af mtools 
package; specifically, the estimation, diagnostic, and forecasting functions. The 
afmtools package would be improved by incorporating other functions related to 
change-point models and tests of unit roots, as well as other important features of 
the models related to long-memory time series. 
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